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The estimation of the ground temperature profile with respect to the depth and time is the key issue in many en- 


gineering applications which use the ground as a source of thermal energy. In the present work, the influence of 


the model components on the calculated ground temperature distribution has been analysed in order to develop an 


accurate and robust model for the prediction of the ground temperature profile. The presented mathematical mod- 


el takes into account all the key phenomena occurring in the soil and on its top surface. The impact of individual 


model elements on the temperature of the soil has been analysed. It has been found that the simplest models and 


the most complex model result in a similar temperature variation over the simulation period, but only at a low 


depth. A detailed analysis shows that a larger depth requires more complex models and the calculation with the 


use of simple models results in an incorrect temperature and a theoretical COP estimation. 
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Introduction 


The recent rapid growth in many engineering and 
geothermal applications which use the ground as a source 
of thermal energy or as a lower heat source has been no- 
ticed. 

These applications typically strongly depend on the 
temperature fluctuations in the near surface region. The 
ground temperature is crucial parameters describing the 
potential of the thermal energy in the soil and directly 
influences the ground base heat pump systems efficiency 
COP (Coefficient of Performance) and SPF (Seasonal 
Performance Factor). An accurate design of the ground 
source system, which uses the ground heat in the winter 
season is a relatively complex task. It requires knowledge 
in the field of heat transport, engineering application, 
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geology as well as environmental protection. 

To provide a sufficient number of data for a proper 
system design, typically, numerical models or simplified 
analytical/empirical models are used [2]. Numerical 
models have the advantage of being able to include com- 
plex physical phenomena, a complicated soil structure or 
realistic boundary conditions. Simulations of this type 
allow a detailed system analysis and a precise selection 
of the system’s element. The designing of a heat ex- 
changer for the lower heat source, without unnecessary 
re-dimensioning, may lead to a significant reduction of 
the investment costs [7,8] however, underestimation, may 
lead to improper system operation, often deep freezing 
and in the worst case, damage the ground heat exchanger. 

In order to estimate precisely the ground potential and 
to conduct a computer analysis, it is necessary to have 
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Nomenclature 

Pas Pw ground, water density, kg/m? 

Cio Cw ground, water specific heat, J/(kg:K) 
k, ground thermal conductivity, W/(m:K) 
t time, s 

Vj the Darcy velocity, m/s 

E emissivity, 


information about the soils geological structure and 
thermo physical properties as well as the accurate boun- 
dary conditions at all surfaces. Without the knowledge of 
those parameters, a computer simulation which includes 
complex models of heat transport in the ground [9] is not 
able to provide correct information for the selection of 
the ground heat exchanger. Also, attention to the initial 
ground temperature profile has to be paid. The analysis 
presented in the literature [10, 11] shows that the soil 
temperature in the near surface region is significantly, 
daily and seasonally, dependent on ambient temperature 
[12], but there are no studies about such models which 
are a combination of different weather elements e.g. solar 
radiation, wind speed or air temperature. 

In the present work, we develop a mathematical and 
numerical model of heat transport in the ground, taking 
into account all the important, or only selected, physical 
phenomena. The analysis of the selected phenomena has 
allowed to decide which heat transfer components are the 
most important for this type of modelling. The estimation 
of the temperature profile [3] and the mean seasonal 
temperature of the ground has allowed to estimate the 
ground thermal potential, which is often one of the most 
important criteria for an economic justification for the 
use of this kind of systems. 


Mathematical Model 


In order to calculate the temperature in the ground, it 
is necessary to construct an accurate mathematical model 
of the occurring transport phenomena [13]. A graphical 
representation of energy balance including the important 
phenomena is given in Figure 1. 

The symbols represent the heat flux exchange by: qconı 
— the natural convection, qcon2 — the forced convection, 
soi — the solar direct and indirect radiation, qier — the 
Earth solar radiation, qeva— the water evaporation from 
the top ground surface. Additionally, q; represents the the 
heat flux related to the groundwater flow, qpa -the heat 
flux related to the phase change and qo - the Earth natural 
heat flux. 

In the present study, the soil consists of several litho- 
logical layers with L; thicknesses and different thermo- 
physical properties (see Table 1). Assuming unsteady 
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o Stefan—Boltzmann constant, W/(m? K$) 
f fraction of evaporation rate 

Th air relative humidity 

Toy reference (sky) temperature, K 

Tair air temperature, K 


phenomena, ground waters flow and phase changes in the 

soil as well as the variations of the soil properties with 

the depth, the equation for the thermal energy transfer 

can be written as follows: 

a C r yT z). 
a (1) 


a 
Cc. — 
Pat a x hy FB 


2. 2. 
k ak a 3 e +S 
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where p,, Py, are the ground and water densities, Cp, 
Cy is the specific heat of the ground and water, k, is the 
thermal conductivity of the ground, t is 
Darcy velocity of the groundwater flow in x direction, v, 
and v, are respectively in y, z direction, and s means the 
source term resulting from the phase change. 


time, v, is the 


Fig. 1 The main elements of the thermal energy balance 


Table 1 Lithological profile of the ground 


Ceiling Floor Thermal conductivity Heat capacity 


[m] [m] [W/mK] [kJ/m*K] 
0 2:3 1.60 2400 
2.3 3.0 1.00 1600 
3.5 4.5 1.20 1700 
4.5 10 1.40 2300 
10 15 0.50 2300 
15 35 2.40 2300 
30 100 2.30 2300 
Mean values 1.485 2128.5 
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It is important to notice that the water content in the 
ground is going to freeze when the temperature is lower 
than 0°C and, in this case, the water property will change 
according to the formula presented in [1]. In the case of 
water freezing the local underground flow velocity 
changes according to the mass conservation equation. 


Elements of the energy balance on the top surface of 
the ground: 


The equation of the energy balance for the ground 
surface depends on the model being considered and it 
may be written, for Model 1, as follows: 

T(x, y,2,t)= f(t) 
and for other models (Model 2-Model 4) as follow: 
- (x Tenen) = fla@l= fla © 
Zz Oz 

The heat flux qw: variable in time, is the difference 
between the heat flux transfer qin into the ground surface 
and the heat flux transfer from the ground surface qour: 

tot = din ~ out (4) 

This study takes into account up to all (in the most 
complex Model 4) the key atmospheric phenomena pre- 
sented in Fig. 1. In accordance with the heat exchange 
methods provided, the heat flux balance for the ground 
top surface can be written as: 


ior = MAX (Geonis Teon2) + Got + Ger +4en ©) 
The first two balance components in eq. (5) represent 
free and forced convection. When the positive tempera- 
ture gradient is detected above the ground surface, the 
conditions for free (natural) convection qconi are met. 
When the wind speed is non-zero forced convection qcon2 
occurs. In the present paper, is was assumed that, at the 
local scale, the earth is a flat horizontal plate and the 
convective heat transport may be determined as: 
dan = Fen Cue Ty) © 
where: Qon-Qcon| OF Acon2 are the heat transfer coeffi- 
cients, in the case of free and forced convection, respec- 
tively and, Tair, is the air temperature. Because of the 
complex model equations for the convective heat transfer 
analysis, usually semi-empirical solutions are used. For 
forced convection, the value of the Nusselt number was 
determined based on the following formula for the lami- 
nar and turbulent flow regime over the flat plate [9]: 
Nu = 0.664Re’?Pr'? (7) 
Nu =0.037Re*°Pr'? for Re>10° (8) 
In the case of natural convection, the following for- 
mula were used [7]: 
Nu =0.54Ra"* for 104<Ra< 10’ (9) 
Nu =0.15Ra'?° for 10’>Ra (10) 
Each time step only one type of convection, the one 
that had higher contributions to the balance (eq. 5) was 
taken into account. 


(2) 
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The second element of the heat balance on the ground 
surface is solar radiation qso; calculated based on the 
measured data. It was assumed in this study that the 
earth's surface can be treated as the approximation of a 
perfect grey body with the emissivity of ¢ = 0.97. 

Due to the differences between ambient temperature 
and ground temperature, the amount of thermal radiation 
heat source qer can be determined from the following 
formula: 

der =€:0-(T(x,y,2,t) -T4,) (11) 

Where Tsw is the reference temperature. The reference 
(or sky) temperature in the present mathematical model 
was calculated with the following formula [12-14]: 

T = 0.0552 T} 


a (12) 
The heat flux related to the evaporation qeva depends 
on several factors (soil moisture, air temperature, veloci- 
ty and humidity) and can be calculated from the follow- 
ing formula [16,17]: 
deo = f Gm [a -T (x, y,0,t)+ b)- r, (a -Ta +b)] 03) 
where f is the fraction of evaporation rate depending 
on the ground moisture and cover, rp is the air relative 


humidity, and a, b are the eq. constants (a=1.73, b=10.23). 


The actual value of air relative humidity was taken from 
a meteorological measurement. The fraction of evapora- 
tion rate f was assumed to be 0.6 which is typical for 
moist soil. 

For the current research, all the models take into con- 
sideration the Earth's heat flux q.=70 mW/m’. 

Based on the heat balance on the surface (eq. 5), four 
different models concerning selected environmental 
phenomena were developed (see Table 2). 


Table 2 Mathematical models for the top surface. 


Assumed conditions on the ground top surface (z=0) 


Model 1 T(x,y,z,t) = KÐ 

Model 2 Tera. flan O] 

Model3 TUK) fg.) +au()+ae(O 

Modig 7H) fg. (0)+ duu (+ dur E) doo O 


dz 
Numerical Model and Results 


In order to solve the presented mathematical model, a 
numerical algorithm was developed in Fortran 90, based 
on the finite volume method. For the convective terms 
and diffusion terms, the second order accurate central 
difference scheme was used CDS). For the spatial deriva- 
tive, three-level second order accuracy methods were 
implemented. The spatial and temporary discretization 
were determined on the basis of the testing calculations, 
which allowed to make the solution independent from 
numerical parameters. In order to speed-up computations 
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the variable time step procedure was used. The geometry 
for analysis including the area with the dimensions of 30 
x 30 x 100 m, in the directions x, y and z, respectively. 
The applied mesh (50 x 50 x 150) was irregular, fine 
close to the surface, where the temperature changes were 
significant and coarse for the depth larger than 30 m. The 
local mean annual temperature Tin; = 9.6°C was used as 
initial and boundary condition for the side geometry 
walls. 

Temporary values of air temperature Tair, wind speed 
magnitude U and solar total radiation lr, were introduced 
into the model from the meteorological measurements for 
the location: Krakow Poland (performed in 2014). In 
order to obtain a statistically steady solution (for a period 
of 1 year), the same weather information was repeated 
for a 10 year calculation period. 
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Fig. 2 Calculated temperature vs. experimental measurement for 
region of Bialystok City [18] and depth z= 1 m (top) 
and z= 2 m (bottom). 


Comparison with experimental measurement: 


The numerical model was validated with the expe- 
rimental measurements of the ground temperature carried 
out in the region of Bialystok, Poland [18]. For the local 
soil thermal properties and the measured air temperature 
the ground temperature at the two different depths (z = 1 
m and z = 2 m) for the time of one year 2001 is presented 
in Figure 3 and compared with the experimental mea- 
surements. 

Taking into account the fact that the exact environ- 
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mental conditions were not known, the air temperature 
was measured at 10 a.m. only (approximately every two 
weeks), while the ground structure except density and 
moisture content were not measured. The agreement be- 
tween the experimental measurements and the numerical 
prediction at two different depths is relatively good. 
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Fig. 3 Comparison of models for a 10-year calculation at the 
depth z= 1 m (top) and z = 5 m (bottom). 


The largest difference occurs at the end of the summer 
season, when very high temperatures occasionally occur 
on the day of the measurement (237th day of the year 
with the temperature of 30.6°C), and this value was used 
until the next measurement was performed. The second 
source of difference between the numerical and the expe- 
rimental results is the initial temperature profile. 

The unsteady solution depends not only on the boun- 
dary conditions but also on the initial conditions. In the 
case of numerical simulations uniform temperature of 
8.5°C was set-up, while in the real soil structure, the 
temperature profile is a result of the historical weather 
phenomena. 

Models analysis: Figure 3 shows the temperature of 
the ground calculated for 10 years with the use of four 
different models (Model 1-Model 4) considered in this 
work. Figure 3 (top) presents the temperature at the depth 
z = ] m. It can be seen that the most complex model 
(Model 4) which is considered here as the reference 
model, could be roughly approximated with all other 
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models. Similar results for the depth z = 5 m, are shown 
in Figure 3(bottom) where the temperature level as well 
as the amplitude of temperature oscillations are much 
lower. Again, also at this depth, all the models follow 
roughly the Model 4. The temperature at this depth is 
going to be statistically steady (for the period of 1 year) 
after approximately 5 years. 

Figure 4 presents temperatures at the four different 
depths z = 0, 1, 2 and 5m and for the last year of simula- 
tions. The calculated top surface temperature depends on 
the mathematical models used. As a consequence, the 
temperatures inside the ground are also slightly different. 
The largest difference at the ground surface (z = 0 m) 
occurs for Model 1. This model is only air temperature 
dependent and “smoothes” the surface temperature. Nei- 
ther the solar radiation nor the wind speed influence the 
surface temperature in this model. It is interesting to no- 
tice that the top surface temperature calculated with 
Modell is always between the results predicted with the 
other models. Additionally, the amplitude of temperature 
fluctuations was always significantly smaller for this 
model. These phenomena are well seen at larger depths. 
The temperature fluctuations predicted with Model 1 are 
smaller also at depths z = 1, 2, 5 m in reference to Model 
4 as well in comparison to the other models. 

Models 2 and 3 result in temperatures which are over- 
estimated in the summer season and in temperatures 
which are very close to the reference models in the win- 
ter season (see Model 3). Finally, Model 3 results in the 
largest temperature fluctuations. 

This model takes into account all the weather pheno- 
mena except for evaporation which significantly stabilis- 
es the temperature fluctuation. All these behaviours are 
very clearly seen at the largest depth z = 5 m, where 
Model 2 predicted a temperature difference mostly in 
reference to Model 4, and the temperature difference is 
almost constant over the whole period and equal to about 
1.5°C. For this depth it can be also clearly seen that the 
amplitude of the temperature fluctuations is the largest 
for Model 3. 

Figure 5 presents the temperature distribution with the 
depth at two different time periods of the year (at the end 
of the year and in the middle of the year). The scope of 
changes regards mainly the area up to ca. 20m of depth. 
In the case of different mathematical models used to cal- 
culate the temperature, the difference occurs also at larg- 
er depths. 

The profiles of temperature for different time periods 
show where the largest temperature variations are ob- 
served. Up to about 1-1.5 m, the freezing zone occurs 
and the ground temperature may drop below 0°C. Ap- 
proximately at this depth area, the largest temperature 
variation is seen. Approximate at the depth z = 1.5 m, the 
first temperature inflection point is seen (see Figure 5 
bottom). The second inflection point is seen at the depth z 
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= 5 m. Depending on the model used for the temperature 
calculation result temperature at different depths and at 
various time period may differ significantly. The largest 
temperature difference between the models was observed 
around the first inflection point in the middle of the year. 
The temperatures predicted with the Model 4 and Model 
2 differ as much as 4°C which is one of the largest dif- 
ferences detected between the models over the whole year. 
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Ground temperature [°C] 


Ground temperature [°C] 


Ground temperature [°C] 


3300 3350 3400 3450 3500 3550 3600 3650 
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Fig. 4 Comparison of mathematical models for last year calcula- 
tion at the depth z =0, 1, 2, 5 m. 
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Fig. 5 The profile of temperature at the end of the simulation 


(top) and in the middle of last year (bottom). 


At the depth of about z = 20 m, the temperature time 
variation is below 0,5°C and this value is influenced only 
by significant weather phenomena variations over a long 
period. 

The mean annual and mean seasonal temperature cal- 
culated for all the models are presented in Table 3 for the 
depth z = 1 m and 5 m. The mean seasonal temperature 
for z = 5 m is in the range from 9.39 to 10.75°C for all 


Table 3 Mean annual and mean seasonal ground temperature. 
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the models and only 9.5 to 9.39°C for Model 1 and Mod- 
el 4. It is worth noticing that the mean values acquired by 
means of the said models overlap highly. The presented 
mean annual and mean seasonal temperature values 
should be treated as the theoretical upper limit value for 
the heat pump or they should be more precise for the 
ground heat exchanger. In the absence of a ground heat 
exchanger (heat pump component) there is no heat trans- 
fer from the ground into the heating system. Typically, 
the heat transfer will decrease the local temperature sig- 


nificantly and influence the overall system performance 
(COP/SPF). 


Conclusions 


In this paper, different mathematical models for 
ground temperature calculations have been studied. The 
influence of the environmental model on the calculated 
ground temperature distribution has been analysed in 
order to develop a fast, accurate and well convergent 
model for the modelling of the ground source heat pump 
system. The presented mathematical model takes into 
account all the key phenomena occurring in the ground 
and on its top surface. The impact of individual compo- 
nents on the temperature of the soil has been analysed. 

Model 4, the most advanced one, taking into account 
all the important phenomena (convection, solar radiation, 
evaporation), has been considered as a reference model 
and in theory, it should match the real conditions in the 
best way. It has been found that none of the models was 
able to reproduce the temperature distribution or temper- 
ature variation with an accuracy similar to the reference 
model. A good approximation with a relatively small 
error has been achieved with the simplest, Model 1. This 
model behaves surprisingly well in comparison to the 
more complex models - Model2 and Model3 — and it is 
computationally cheap and numerically stable. 


Model 1 Model 2 Model 3 Model 4 
mean annual vs conduction+ 
conduction conduction + . 
mean seasonal convection+ 
conduction + convection+ S 
(1 Sep. — 31 Mar.) i A thermal radiation + 
convection thermal radiation : 
evaporation 
T [°C] 9.78 9.64 9.64 9,64 
z=0m a 
T [°C] 16.34 16.23 16.23 16.23 
i T [°C] 9.75 11.36 10,62 9.61 
z=Ilm 2 
T [°C] 14.30 15.98 16.85 15.11 
T [°C] 9.74 11.27 10,57 9.6 
z=2m _ 
T [°C] 12.56 14.13 14.48 13.03 
T [°C] 9.71 10.98 10,39 9.59 
z=5m B 
T [°C] 9.5 10.75 10.21 9.39 
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It has been found that the simplest models and the 
most complex model result in different temperature fluc- 
tuations over the simulation period, but the mean value 
may not differ significantly. A detailed analysis shows 
that a larger depth typically required more complex mod- 
els and a calculation with the use of simple models re- 
sults in an incorrect maximum and minimum as well as 
mean temperature and in consequence theoretical COP, 
SPF estimation. 

For Model 1, the mean seasonal temperature predic- 
tion at the surface z = 0 m differs only by 0.14°C (see 
Table 3) in reference to Model 4. This means that, under 
the conditions taken into account in this work, the air 
temperature and the wind speed are good representative 
factors for all the environmental phenomena, regardless 
of the type of ground (layer thickness, lithology or satu- 
ration). 

In any case, the temperature for the depth z>30m is 
almost model independent and, with a good approxima- 
tion, equals mean annual air temperature. In order to 
obtain a statistically steady temperature at larger depths, 
a much longer computational time (about 50 years) is 
required. 
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